home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / SPARSE.DEM < prev    next >
Text File  |  1991-04-29  |  2KB  |  74 lines
  1. PROGRAM d2r10(input,output);
  2. (* driver for routine SPARSE *)
  3. CONST
  4.    n=20;
  5. TYPE
  6.    glnarray = ARRAY [1..n] OF real;
  7. VAR
  8.    i,ii : integer;
  9.    rsq : real;
  10.    b,bcmp,x : glnarray;
  11.  
  12. PROCEDURE asub(xin: glnarray; VAR xout: glnarray; n: integer);
  13. (* Programs using routine ASUB must define the type
  14. TYPE
  15.    glnarray = ARRAY [1..n] OF real;
  16. in the main routine. *)
  17. VAR
  18.    i: integer;
  19. BEGIN
  20.    xout[1] := xin[1]+2.0*xin[2];
  21.    xout[n] := -2.0*xin[n-1]+xin[n];
  22.    FOR i := 2 to n-1 DO BEGIN
  23.       xout[i] := -2.0*xin[i-1]+xin[i]+2.0*xin[i+1]
  24.    END
  25. END;
  26.  
  27. PROCEDURE atsub(xin: glnarray; VAR xout: glnarray; n: integer);
  28. (* Programs using routine ATSUB must define the type
  29. TYPE
  30.    glnarray = ARRAY [1..n] OF real;
  31. in the main routine. *)
  32. VAR
  33.    i: integer;
  34. BEGIN
  35.    xout[1] := xin[1]-2.0*xin[2];
  36.    xout[n] := 2.0*xin[n-1]+xin[n];
  37.    FOR i := 2 to n-1 DO BEGIN
  38.       xout[i] := 2.0*xin[i-1]+xin[i]-2.0*xin[i+1]
  39.    END
  40. END;
  41.  
  42. (*$I MODFILE.PAS *)
  43. (*$I SPARSE.PAS *)
  44.  
  45. BEGIN
  46.    FOR i := 1 to n DO BEGIN
  47.       x[i] := 0.0;
  48.       b[i] := 1.0
  49.    END;
  50.    b[1] := 3.0;
  51.    b[n] := -1.0;
  52.    sparse(b,n,x,rsq);
  53.    writeln('sum-squared residual:',rsq:15);
  54.    writeln;
  55.    writeln('solution vector:');
  56.    FOR ii := 1 to (n DIV 5) DO BEGIN
  57.       FOR i := (5*(ii-1)+1) to (5*ii) DO write(x[i]:12:6);
  58.       writeln
  59.    END;
  60.    IF ((n MOD 5) > 0) THEN BEGIN
  61.       FOR i := 1 to (n MOD 5) DO write(x[5*(n DIV 5)+i]:12:6)
  62.    END;
  63.    writeln;
  64.    asub(x,bcmp,n);
  65.    writeln;
  66.    writeln('press RETURN to continue...');
  67.    readln;
  68.    writeln('test of solution vector:');
  69.    writeln('a*x':9,'b':12);
  70.    FOR i := 1 to n DO BEGIN
  71.       writeln(bcmp[i]:12:6,b[i]:12:6)
  72.    END
  73. END.
  74.